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Abstract. System identification or mathematical modelling is utilised in the aerospace com- 
munity for the development of simulation models for robust control law design. These models 
are often described as linear, time-invariant processes and assumed to be uniform throughout 
the flight envelope. Nevertheless, it is well known that the underlying process is inherently 
nonlinear. The reason for utilising a linear approach has been due to the lack of a proper set of 
tools for the identification of nonlinear systems. Over the past several decades the controls and 
biomedical communities have made great advances in developing tools for the identification of 
nonlinear systems. These approaches are robust and readily applicable to aerospace systems. 

In this paper, we show the application of one such nonlinear system identification technique, 

TM 

structure detection, for the analysis of F-15B QuietSpike aeroservoelastic flight test data. 
Structure detection is concerned with the selection of a subset of candidate terms that best 
describe the observed output. This is a necessary procedure to compute an efficient system 
description which may afford greater insight into the functionality of the system or a simpler 
controller design. Structure computation as a tool for black-box modelling may be of critical 
importance for the development of robust, parsimonious models for the flight-test community. 
Moreover, this approach may lead to efficient strategies for rapid envelope expansion which 
may save significant development time and costs. 

TM 

The objectives of this study are to demonstrate via analysis of F-15B QuietSpike aeroser- 
voelastic flight test data for several flight conditions (Mach number) that (i) linear models are 
inefficient for modelling aeroservoelastic data, (ii) nonlinear identification provides a parsimo- 
nious model description whilst providing a high percent fit for cross-validated data and (iii) the 
model structure and parameters vary as the flight condition is altered. 

*This work was prepared as part of the author’s official duties as an employee of the U. S. Government and in 
accordance with 17 U.S.C. 105, is not available for copyright protection in the United States. NASA is the owner 
of any foreign copyright that can be asserted for the work. Copyright©2007 by NASA. 
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1 Introduction 


System identification, or black-box modelling, is a critical step in aircraft development, analy- 
sis and validation for flight worthiness. The development and testing of aircraft typically takes 
many years and requires considerable expenditure of limited resources. One reason for lengthy 
development time/costs is assuming the underlying system is linear and invariant throughout 
the flight envelope. This assumption is related to having an inadequate knowledge of an ap- 
propriate model type or structure to use for parameter estimation. Selection of an insufficient 
model structure may lead to difficulties in parameter estimation, giving estimates with signif- 
icant biases and/or large variances [1]. This often complicates control synthesis or renders it 
infeasible. The power of using structure detection techniques as a tool for model development 
(i.e. black-box modelling) is that it can provide a parsimonious system description which can 
describe complex aeroelastic behaviour over a large operating range. Consequently, this pro- 
vides models that can be more robust and, therefore, reduce development time. 

Moreover, when studying aeroelastic systems it may not be practical to assume that the exact 
model structure is well known a priori. In aerospace systems analysis one of the main ob- 
jectives is not only to estimate system parameters but to gain insight into the structure of the 
underlying system. Therefore, structure computation is of significant relevance and importance 
to modelling & design of aircraft and aerospace vehicles. Structure computation may indicate 
deficiencies in an analytical model and could lead to improved modelling strategies and also 
provide a parsimonious, black-box, system description for control synthesis [2]. 

For linear systems modelling a commonly used approach for determining model structure is the 
minimum description length (MDL) proposed by Rissanen [3]. MDL was specifically devel- 
oped to overcome some of the inconsistencies of the Akaike’s information criterion (AIC) [4], 
i.e. its variance does not tend to zero for larger sample sizes, N. 

Recently, the bootstrap method has been shown to be a useful tool for structure detection of non- 
linear models [5, 6, 7]. The bootstrap is a numerical method for estimating parameter statistics 
which requires few assumptions [8]. The conditions needed to apply bootstrap to least-squares 
estimation are quite mild; namely, that the errors be independent, identically distributed, and 
have zero-mean. 

Identification of flight dynamics are not well understood in the subsonic, transonic and su- 
personic regimes. Therefore, in this paper, we investigate whether (i) a linear or nonlinear 
model best represents the observed data and (ii) the system is invariant during envelope expan- 

TM 

sion (varying Mach number). The data analysed in this study is from the F-15B QuietSpike 
flight test program which was a collaborative effort between Gulfstream Aerospace Corpora- 
tion (Savannah, Georgia) and NASA Dryden Flight Research Center [9, 10, 11]. The results 
show that (i) linear models are inefficient for modelling aeroservoelastic data, (ii) nonlinear 
identification provides a parsimonious model description whilst providing a high percent fit for 
cross-validated data and (iii) the model structure and parameters vary as the flight condition is 
altered. 
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2 NARMAX Model Form 


The dynamic behavior of many nonlinear systems can be represented as a discrete polynomial 
which expands the present output value in terms of present and past values of the input signal 
and past values of the output signal [12, 13, 14]. A system modelled in this form is popularly 
known as a NARMAX (Nonlinear AutoRegressive, Moving Average exogenous) model and is 
linear- in-the -parameters. 

Recently, Kukreja el a!. [15] showed that NARMAX identification is well suited to describing 
aeroelastic phenomena. The NARMAX structure is a general parametric form for modeling 
nonlinear systems [12]. This structure describes both the stochastic and deterministic compo- 
nents of nonlinear systems. Many nonlinear systems are a special case of the general NARMAX 
structure [16]. In this paper, we focus on a special class of NARMAX models; nonlinear poly- 
nomial models. The polynomial NARMAX structure models the input-output relationship as a 
nonlinear difference equation of the form 

z(n) = f l [ z ( n ~ I),-" , z(n -n y ),u(n),- ■ ■ , u(n - n u ) , e(n - l),--- , e(n - n e )\ +e(ri). (1) 

/ denotes a nonlinear mapping, l is the order of the nonlinearity, u is the controlled or exogenous 
input, z is the measured output, and e is the uncontrolled input or innovation. This nonlinear 
mapping may include a variety of nonlinear terms, such as terms raised to an integer power, 
products of present and past inputs, past outputs, or cross-terms. In general, the nonlinear 
mapping, /, can be described by a wide variety of nonlinear functions such as sigmoids or 
splines [17, 16]. This system description encompasses many forms of linear and nonlinear 
difference equations that are linear-in-the -parameters. 

Identifying a NARMAX model requires two things: (1) structure detection and (2) parameter 
estimation. Structure detection can be divided into: (la) model order selection and (lb) se- 
lecting which parameters to include in the model. We consider model order selection as part 
of structure detection since, theoretically, there are an infinite number of candidate terms that 
could be considered initially. Establishing the model order, then, limits the choice of terms to 
be considered. For the NARMAX model, the system order is defined to be an ordered tuple as 

O = [n u n z n e 1} (2) 

where n u is the maximum lag on the input, n z the maximum lag on the output, n e the maxi- 
mum lag on the error and l is the maximum nonlinearity order. Note that for non-polynomial 
NARMAX models, l may be simply replaced by a nonlinear mapping of some specified class. 
In this paper, we assume that the system order is known. 


3 Structure Detection 


The structure detection problem is that of selecting the subset of candidate terms that best de- 
scribes the output. Therefore, the parametrisation of a system is still further reduced by deter- 
mining which of the components are required. The maximum number of terms in a NARMAX 
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model with n z , n u and n e dynamic terms and Zth order nonlinearity is 


i 

p = + 1 \Pi 

2=1 


Pi-ii;n u + n z + n e + i) 
: ,Po = 1- 

i 


( 3 ) 


As a result, the number of candidate terms becomes very large for even moderately complex 
models making structure detection difficult. We define the maximum number of terms, p, as 
the number of candidate terms to be initially considered for identification. Due to the excessive 
parameterisation (the curse of dimensionality), the structure detection problem often leads to 
computationally intractable combinatorial optimisation problems. 


4 Time Series 


The data considered in this paper is time series since the input signal, u(n), is assumed to be 
zero or constant. Time series analysis is used when inputs are not available to the experimenter, 
or where it is unclear which signals are inputs and which are outputs [18]. Models arising from 
times series data can have many forms [1]. However, in our treatment of the data the ARMA 
and NARMA model class are of practical significance. 

This special case of the general NARMAX model (Eqn. 1) can be written as 

z(n) = f[ z ( n - I),"’ ,z(n-n z )--- , e(n - n e )] + e(n) (4) 

where we redefine the model order for this model set as 

0 = [n z n e l\. (5) 

The maximum number of candidate terms in a model (Eqn. 4) with n z and n e dynamic terms 
and Zth order nonlinearity is 

i 

P = 5^P* + 1; (6) 

2=1 

Pi-i(n z +n e + i) 

Pi = : , Po = I- 

I 


Note that ARMA models can be estimated using the Yule-Walker equations or the instrumental 
variable (IV) estimator to avoid estimating the MA part [1], This is the approach taken in this 
paper. For a NARMA model the NMA part must be modelled. For nonlinear systems, output 
additive noise can produce multiplicative terms between input, output and itself. To compute 
unbiased parameter estimates a noise model (i.e.NMA) needs to be estimated [19]. 
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5 Structure Detection Methods 


With the model types defined for the flight test data available for analysis, we describe two 
approaches applicable to these model classes. The first is appropriate for AR models whilst the 
second for NARMA models. 


5.1 Minimum Description Length 

A commonly used technique in linear system identification to determine model structure is 
MDL [3]. MDL was specifically developed to overcome some of the inconsistencies of the 
Akaike’s information criterion (AIC) [4], i.e. its variance does not tend to zero for larger sample 
sizes, N. 

The number of parameters necessary to reproduce an observed sequence {zi , . . . , zn} of a time 
series depends on the model and parameters assumed to have generated the data [3]. The MDL 
technique finds the model which minimizes the description length and thereby computes an 
estimate of model order [3]. 

Binary prefix codes are used to encode data strings. These data strings can be made up of 
symbols, parameters, numbers, etc. It is known that the average length of a code word is 
bounded by Shannon’s theorem [3]. Therefore, it is possible to write [3] 

p(x)L(x ) > — ^^p( x ) logp(x) (7) 

X X 


where L(x) is the length of the code word (i.e. length of parameter vector 6) and p(x) is the 
probability of x. It is also possible to write 

L(z\x, 6) = — \ogp(z\x, 9) (8) 

where L(z \x, 9) is known as the log-likelihood function (to be maximised). Let 0 denote the 
value of the parameter which maximises the likelihood and thus minimises the parameter vector 
length (i.e. code word length) L(y\x,9). Since 9 can only be encoded up to a certain precision, 
the code word length, L(y\x,9), becomes longer than the desired minimum L(z\x. 9). given 
noise considerations. Let the precision be 5 = 2~ q where q is the number of bits used for 
encoding the parameter. It is possible to save on the code word length if q is small. However, 
the result is a loss in precision. The optimal precision depends on the size of the observed data 
via — log 5 = 0.5 log N, and hence the total code word length for k parameters is given by the 
MDL, 


MDL(fc) = — log [maximised likelihood] + -A; log A" 
which, for an AR(n 2 ) model gives 

77 / 

MDL(n^) = log [maximised likelihood] + — log N. 


(9) 


( 10 ) 
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5.2 Bootstrap 


Recently, the bootstrap has been shown to be a useful tool for structure detection of nonlin- 
ear models [6]. The bootstrap is a numerical method for estimating parameter statistics which 
requires few assumptions [8]. The conditions needed to apply bootstrap to least-squares esti- 
mation are quite mild; namely, that the errors be independent, identically distributed, and have 
zero-mean. 

The bootstrap is a technique to randomly reassign observations which enables re-estimates of 
parameters to be computed. This randomisation and computation of parameters is done numer- 
ous times and treated as repeated experiments. In essence, the bootstrap simulates a Monte- 
Carlo analysis. For structure computation, the bootstrap method is used to detect spurious 
parameters, those parameters whose estimated values cannot be distinguished from zero. 

Application of an appropriate t 2 estimator to measured data gives the model response Z, residu- 
als £ and parameter estimate 6. The bootstrap assumes that these residuals, £ — [ei, £ 2 , • • • , £jv], 
arise from an unknown distribution, V. By randomly resampling these residuals, with replace- 
ment, it is possible to generate aresampled version of the prediction errors, £* = [£*, i 2 , ■ ■ ■ , £^], 
whose distribution estimates V. The resampling procedure for each £* involves randomly se- 
lecting from £ with an equal probability associated with each of the N elements. For example, 
a possible resampled version of the errors for N = 5 is £* = [f 4 , £\, £4, £ 2 , £3]. The star notation 
indicates that £* is not the calculated error £ but rather a resampled version of it. These resam- 
pled errors are added to the model response to generate a bootstrap replication of the original 
data, 


z* = ^ > Z£ e + £*. 


(ii) 


A new bootstrapped parameter estimate 6 is obtained from this bootstrapped data Z*. This 
procedure is repeated B times to provide a set of parameter estimates from the B bootstrap 
replications, 


©* 



( 12 ) 


Parameter statistics can then be easily computed from ©* by forming percentile intervals at a 
chosen level of significance, a. 


Structure detection can provide useful process insights that can be used in subsequent develop- 
ment or refinement of physical models. Therefore, in the sequel, we investigate the applicabil- 
ity of MDL and the bootstrap to experimental aircraft data. Specifically, MDL and bootstrap 
methods are used as structure detection tools to assess whether the (i) underlying data is best de- 
scribed by a linear time-invariant (LTI) or nonlinear model and (ii) model structure is invariant 
during envelope expansion. 


6 Experimental F-15B QuietSpike' 1 Data 


The MDL and bootstrap technique were applied to experimental flight test data from the F-15B 

TM 

QuietSpike project by Gulfstream Aerospace Corporation (Savannah, Georgia) and NASA 
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Dryden Flight Research Center [9, 10, 11]. The data analysed for this study used structural 

TM 

accelerometer response output of the QuietSpike boom tip when fully extended. 


6.1 Procedures 


Flight data was gathered during subsonic, transonic and supersonic flutter clearance of the F- 

TM 

15B QuietSpike (see Fig 1). This paper considers accelerometer data measured during pilot 



Figure 1: Flight test article in extended configuration, 
induced pitch-raps at Mach 0.85, 0.95 and 1.40 at 12,192 m (40, 000 ft). The output was taken 

TM 

as the response of an accelerometer mounted near the QuietSpike boom tip (see Fig. 2). Data 
was sampled at 400 Hz. For analysis, the recorded flight test data was decimated by a factor of 
2, resulting in a final sampling rate of 200 Hz. 

For identification of a linear model an arbitrarily large AR model of fiftieth order (n z = 50) was 
posed for estimation and the MDL technique used to determine the optimal linear model. For 
identification of a nonlinear model a model order with fourth-order output and error dynamics 
and second-order nonlinearity, O = [4,4,2], was used. A model with fourth-order dynamics 
was selected because it has been observed that aeroservoelastic structures are well defined by 
a fourth-order LTI system [20]. The nonlinearity order was chosen as second-order because 
empirical results showed models of higher nonlinear order were not efficient to describe the 
data sets available for analysis. This gave a full model description with 45 candidate terms. 
The nonlinear model was identified applying the bootstrap approach. For the bootstrap method, 
B = 100 bootstrap replications were generated to assess the distribution of each parameter. For 
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Primary 


Backup 



TM 

Figure 2: QuietSpike sensor locations. 



all three techniques, each parameter was tested for significance at the 95% confidence level. 

For both linear and nonlinear identification, Table 1 shows the number of data points available 
for each flight condition. The estimation data sets were from accelerometer response measure- 


Altitude 12,192 m (40, 000 ft) 


Mach Number 

0.85 

0.95 

1.40 

Estimation Data 

572 

400 

504 

Validation Data 

572 

400 

504 


Table 1: Data points available at each flight condition. 

ments of the primary sensor on the boom tip and the cross-validation data sets were of the 
backup sensor at the same location. 


6.2 Results 


TM 

The results of identifying the F-15B QuietSpike data are presented. Fig. 3 shows the output 
data sets used for this analysis. The data represents structural accelerometer response (primary 
sensor boom tip) used to compute the system structure. 
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Structural Accelerometer Response Structural Accelerometer Response 




(a) Mach 0.85 (b) Mach 0.95 


Structural Accelerometer Response 



(c) Mach 1.40 


Figure 3: Estimation data: Recorded structural accelerometer response to stick raps. 
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Eqns. 13-15 depict the model structure computed by the bootstrap method 

Mach 0.85, Altitude 40,000 ft. (12,192 m) (13) 

z(n) = 6\z{n — 1) + fa z(n — 2) + 9 3 z(n — 3) + faz 2 {;n — 4) + fafan — 1) 

+ fafan — 2) + 9 7 e(n — 3) + faz{n — 4 )e(n — 4) + 9g£ 2 {n — 4) + fan) 

Mach 0.95, Altitude 40,000 ft. (12,192 m) (14) 

z(n) = diz(n — 1) + d 2 z(n — 2) + d 3 z(n — 3) + faz(n — 3)z(n — 4) 

+ 9 5 z 2 (n — 4) + i? 6 e(n — 1) + S 7 s(n — 2) + fae{n — 3) 

+ '9 9 z(n — 3 )e(n — 4) + d w z(n — 4 )fan — 3) + 9 n £(n — 3 )e(n — 4) 

+ #12 z(n — 4 )fan — 4) + #1 3 £ 2 (n — 4) + e(n) 

Mach 1.40, Altitude 40,000 ft. (12,192 m) (15) 

z(n) = faz(n — 1) + fa z(n — 2) + fa z(n — 3) + fa z(n — 1 )z(n — 4) 

+ fafan - 1) + fa fan - 2) + fafan - 3) + faz{n - 1 )e(n - 4) 

+ faz(n — 4)e(n — 1) + fao£(n — l)e(n — 4) + e(n). 

Eqns. 13 - 15 represent the computed model structure for flight conditions Mach 0.85, 0.95 and 
1.40, respectively. The computed model structures are represented by a combination of linear 
and nonlinear, lagged input-output terms and contain 9, 13 and 10 terms for Mach 0.85, 0.95 
and 1.40, respectively. Hence, the bootstrap technique successfully produced a parsimonious 
model description from the full set of 45 candidate terms. 

For AR (linear) model identification using MDL to compute structure, the estimated models 
were of order n z = 42, 44, and 46 for Mach 0.85, 0.95 and 1.40, respectively. These models 
are not shown since they are simply a dynamic expansion of the output up to the order stated. 
However, for cross-validation data, we show the model fit of these linear models compared to 
the cross-validation fit obtained with the NARMA models (see Fig. 4). 

Fig. 4 shows the predicted output for the cross-validation data sets for the identified structures 
((a): Eqn. 13 and AR (n z = 42), (b): Eqn. 14 and AR(>i 2 = 44) and (c): Eqn. 15 and AR(n z = 
46)). Each panel displays the full time history of the predicted output of the linear and nonlinear 
models superimposed on top of the measured output. For Mach 0.85 (Fig. 4(a)) the predicted 
output of the linear and nonlinear models account for over 91% and 95% of the measured 
outputs variance, respectively. For Mach 0.95 (Fig. 4(b)) the predicted output of the linear and 
nonlinear models account for over 92% and 97% of the measured outputs variance, respectively. 
For Mach 1.40 (Fig. 4(c)) the predicted output of the linear and nonlinear models account for 
over 91% and 96% of the measured outputs variance, respectively. 

The results demonstrate that although the AR models contain many more terms to explain the 
underlying process they still offer a lower percent fit compared to the nonlinear model at the 
cost of model complexity (higher order) which often leads to more complex control synthesis. 
The nonlinear models contain only a few terms and were capable of explaining a larger per- 
cent of the output variance. For these data sets the results show linear models are inefficient 
for accurate modelling of aeroservoelastic data. These results show a nonlinear identification 
approach offers a parsimonious system description whilst providing a high percent fit for cross- 
validated data. Moreover, the results illustrate the need to vary model structure for different 
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Cross-Validated Structural Accelerometer Response 



Cross-Validated Structural Accelerometer Response 



(a) Mach 0.85 


(b) Mach 0.95 


Cross-Validated Structural Accelerometer Response 



(c) Mach 1.40 


Figure 4: Cross-validation data: Predicted linear and nonlinear model accelerometer response 
of z-tip longitudinal sensor superimposed on top of measured accelerometer output. 
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flight condition. 


7 Discussion 


Experimental results demonstrate that structure computation as a tool for black-box modelling 
may be useful for the analysis of dynamic aircraft data. The bootstrap successfully reduced the 
number of regressors posed to aircraft aeroelastic data yielding a parsimonious model structure 
for each data set. Additionally, these parsimonious structures were capable of predicting a large 
portion of the cross-validation data, collected on a backup sensor at a similar location. However, 
for linear analysis, the MDL approach was not able to reduce the model order (structure) as 
well and yielded a more complex system description. Whilst these linear models have higher 
complexity (degrees of freedom), they provided a model predictive capability that explained a 
smaller percent of the observed output variance. This find indicates a linear model may not be 
appropriate to describe aeroservoelastic data. A higher percent fit offered by the parsimonious 
nonlinear models suggests that the identified structures and parameters explain the data well. 
Using percent fit alone as an indicator of model goodness could lead to incorrect interpretations 
of model validity. Nevertheless, in many cases for nonlinear models this may be the only 
indicator that is readily available. 

In this work, the results show that whilst the linear dynamics remained invariant for all flight 
conditions available for analysis, the nonlinear dynamics changed as Mach number increased. 
For Mach 0.85 the model (Eqn. 13) displayed a mildly nonlinear process which physically 
makes sense since the airflow is mainly subsonic. When the Mach number was increased to 
0.95 the model (Eqn. 14) demonstrated a richer nonlinear dynamic description which is likely 
due to embedded shock formations in the transonic regime. For Mach 1.40 the model (Eqn. 15) 
displayed a mildly nonlinear process again which physically makes sense since in this regime 
the shocks become fixed. It is difficult to make definitive comments on the underlying physics 
responsible for this behaviour without extensive analysis of different flight conditions. The 
important points to note are, this study suggests (i) nonlinear models are appropriate to describe 
the dynamics behaviour of advanced aircraft and (ii) models describing aircraft dynamics vary 
with flight condition. This suggests nonlinear modelling may afford a robust and parsimonious 
system description over a larger operating regime and models used for prediction (e.g. control) 
should not be invariant for all flight conditions. This may hold significant implications for 
aircraft development. 

For this study, only a polynomial mapping with fourth-order output and error lag was used as 

TM 

a basis function to explain the nonlinear behaviour of the F-15B QuietSpike data. Clearly, 
different basis functions and a higher dynamic order (lag-order) should be investigated to deter- 
mine if another basis can produce accurate model predictions with reduced complexity. More- 
over, further studies are necessary to evaluate whether the model structure is invariant under 
different operating conditions, such as altitude, and model parameterisations. 

This study illustrates the usefulness of structure detection as an approach to compute a par- 
simonious model of a highly complex nonlinear process, as demonstrated with experimental 
data of aircraft aeroelastic dynamics. Moreover, analysis of flight test data can provide useful 
process insights that can be used in subsequent development or refinement of physical models. 
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In particular, morphological models are based on assumptions (e.g. these effects are important 
and those are negligible) which may be incorrect [21, 22]. A structure computation approach to 
model identification may help uncover such surprises. 


8 Conclusions 


Results show that linear models are inefficient for modelling aeroservoelastic data and nonlin- 
ear identification provides a parsimonious model description whilst providing a high percent fit 
for cross-validated data. Moreover, the results demonstrate that model structure and parameters 
vary as the flight condition varies. These results may have practical significance in the analysis 
of aircraft dynamics during envelope expansion and could lead to more efficient control strate- 
gies. In addition, this technique could allow greater insight into the functionality of various 
systems dynamics, by providing a quantitative model which is easily interpretable. 
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